c	muladd(a,b,c,m)
c	int a,b,c,m[2];
c	computes m(1) = (a*b + c) mod 2^30 and m(2) = (a*b + c)/2^30

	subroutine muladd(a,b,c,m)
	integer a,b,c,m(2)
	double precision frad1
	frad1 = 1.0d0/1073741824.0d0
	m(1) = and((c + a*b) , 1073741823)
	m(2) = ifix(0.25+frad1*( (dble(c-m(1)))+(dble(a))*(dble(b)) ))
	return
	end

c	divrem(a,b,c,m)
c	int a,b,c,m[2];
c	computes m[1] = a*2^30+b mod c and m[0] = (a*2^30+b)/c
c	assumes a,b,c are all unsigned and less than 2^30
c	uses double length registers

	subroutine divrem(a,b,c,q)
	integer a,b,c,q(2)
	integer lr, lq
	double precision frad
	integer m(2), lo
	lo = 0
	frad = 1073741824.0d0
	lq = ifix((((frad*dble(a)) + dble(b)) / dble(c)))
	call muladd(lq,c,lo,m)
	lr = 1073741824*(a - m(2)) + (b - m(1))

c	if (lr < 0) { do {lq--;} while ((lr += c) < 0);}
c	else	    { while (lr >= c) {lr -= c; lq++;};}
	if (lr .ge. 0) goto 2000
200	lq = lq - 1
	lr = lr + c
	if (lr .lt. 0) goto 200
2000	if (lr .lt. c) goto 3000
	lr = lr - c
	lq = lq + 1
	goto 2000
3000	continue
	q(1) = lq
	q(2) = lr
	return
	end

c	mod_mult(a,b,c)
c	int a,b,c;
c	returns a*b % c
c	assumes a,b,c are less than 2^30

	integer function mod_mult(a,b,c)
	integer a,b,c
	integer m(2), q(2), lo
	lo = 0
	call muladd(a,b,lo,m)
	call divrem(m(2),m(1),c,q)
	mod_mult = q(2)
	return
	end
